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Abstract 


In ultrasonic non-destructive testing and evaluation, the different features of ultrasonic 
waves, namely frequency, amplitude, time of flight etc can be used for detection and 
sizing of internal defects. Ultrasonic Tomography is used to characterize a defect or 
inhomogeneity by time of flight information. The accuracy of reconstruction depends on 
the accuracy of time of flight estimation. The present work is concerned with 
tomographic reconstruction of square, circular and elliptical defects using simulated time 
of flight data. Three combinations of base material and insert material, namely, (1) Steel 
with Aluminum insert, (2) Steel with Nickel insert and (3) Steel with Brass insert have 
been taken in this present work to study the quality of reconstruction with the relative 
variation in longitudinal sonic velocity between the base and the insert material. SnelTs 
law has been used to ascertain the path taken up by the ultrasonic ray in the presence of a 
refraction artifact. Standard edge detection techniques have been applied to find out the 
interface for refraction in the presence of a curved defect. From the nature of the 
reconstructed images the optimum range for the relaxation parameter has been identified 
within which there is a good balance between the quality of reconstruction and also the 
computer run time. All the reconstructions have been obtained using Mayinger ART 
algorithm. The effect of relaxation parameter on the rate of convergence of the iterative 
algorithm has also been studied. Tables providing values of normalized RMS errors have 
been given for all the combinations of the base and insert materials and for all defect 
geometries. The reconstructed images obtained using curved ray approach have been 
compared with the same obtained using straight ray assumption. Another table furnishing 
the values of the normalized RMS errors in the case of straight ray approach has been 
given. Finally, a table showing a comparison of the times taken to execute the code using 
curved ray approach and without using curved ray approach has been provided. 
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Chapter 1 


Introduction 

1.1 Introduction 

Non-destructive testing and evaluation methods play a pivotal role in quality assurance of 
structural members during their manufacturing stage and later during their operational 
period. Being one of the most commonly used non-destructive testing methods, ultrasonic 
testing is rapidly developing in the recent years. This method uses ultrasonic "waves for 
detection and sizing of internal flaws present in materials. In this method, piezo-electric 
transducer (probe) generates ultrasonic waves, which propagate in the elastic medium and 
are detected either by the same or by a different transducer. 

In general, the transit time of the acoustic echo wave is used to measure the 
distance of a flaw or insertion firom the probe and the amplitude pf the reflected wave is 
used to size the flaw. In the pulse-echo technique, a single probe is used to transmit and 

I 

receive the signal. Another technique, called through transmission technique, employs 
two separate transducers (a transmitter and a receiver) to transmit and receive the wave. 
This technique requires access to both sides of the specimen, as the probes are placed on 
opposite sides of the specimen. 

Ultrasonic signals can be displayed in A-scan, B-scan or C-scan formats. An A- 
scan indicates a variation of the signal amplitude with time at a point. The familiar screen 
image is an A-scan presentation, which furnishes one-dimensional description of a given 
test point. In the case of B-scan, the displayed signal is the trace of time versus linear 
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position. In the case of two-dimensional scanning of a test-piece, the results can be 
represented by means of a C-scan. This furnishes a top-view of the scanned surface. C- 
scan systems are especially computer-controlled. Ultrasonic systems displaying both B- 
scan and C-scan formats are also called P-scan systems. 

Tomography is a procedure where internal body images at a predetermined plane 
are recorded by means of the tomograph, a computer-driven machine that builds the 
image from multiple x-ray measurements. It is also called body section roentgenography. 
It can produce sophisticated cross-sectional or three-dimensional pictures of the internal 
organs, which can be used to identify abnormalities, for example, tumors. There are 
numerous applications for reconstructing images from their projections. The spatial 
resolution and tissue discrimination capabilities of this technique far surpass what could 
be achieved with conventional x-ray imaging. Image reconstruction from projections is 
also used in nuclear medicine to map the distribution of the concentration of gamma-ray 
emitting radionuclides in a given cross section of the human body. Another medical 
application area is computed tomographic imaging with ultrasound. Due to the radiation 
hazards associated with x-rays, ultrasonic imaging is gaining considerable attention. Non- 
medical applications in which images can be reconstructed from projections include 
radioastronomy, optical interferometry, electron microscopy and geophysical exploration. 

In contrast with X-rays, tomographic imaging with ultrasound is made difficult 
due to the fact that the rays of sound energy do not necessarily travel in a rectilinear 
fashion when an interface is encountered. This phenomenon of refraction while using 
ultrasound as a source has posed major challenges in the field of tomographic 
reconstraction. To account for refraction, the acoustic analogue of Snell’s law is widely 
used. 

This study is concerned with tomographic reconstruction of isotropic specimens 
with artificially implanted defects. For the above-mentioned reconsfruction, the 
propagation characteristics (speed or amplitude) of the acoustic wave are measured. 
These propagation characteristics, namely, speed and amplitude are directly related to the 
material property of interest. Therefore, the properties, which can be measured 
experimentally, can often be approximated as line integrals along the path of the 
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propagating disturbance. For acoustic experiments in isotropic media, one possible 

t 

formulation of the problem is in terms of transit time between the source and the receiver. 

1.2 Literature Survey 

Theoretical studies on the elastic wave propagation found in literature are manifold, most 
of which deal with isotropic media, only a few are concerned with wave propagation in 
anisotropic media. Extending the scope of interest to the wave scattering problem in 
anisotropic materials for the' purpose of location and classification of defects and other 
inhomogeneities, respective studies are approximate in character and mostly confined to 
particular geometries and material symmetry. Chen [1] first presented a coordinate free 
approach, which is useful in dealing with lay structures. 

Synge [2] developed the theory of slowness surfaces. After that Musgrave [3] 
introduced the idea of superimposing the slowness surface on the geometric space, which 
contains the incident and emerging waves. 

Previously, modelling techniques for inhomogeneous materials have been 
developed with emphasis on inspection on austenitic welds and claddings. Although the 
numerical methods such as finite difference, finite element or finite integration provide 
good results, they suffer from the inherent requirement of large computation times and 
hardware requirements 

Leander [4] developed the relation between the wavefront speed, which is the 
speed of propagation of a front of a pulse. Crandell’s [5] work was based on the use of 
slowness surfaces to represent wave reflections. 

Very few quantitative results have been published on the reflection and 
refraction of plane elastic waves from an interface in a general anisotropic media. 
Musgrave and Fedorov [6] summarized the general theoretical results for the reflection- 
refraction problems in their comprehensive books. Auld [7] discussed several examples 
for planes of symmetry. Henneke [8] and McNiven et al [9] reviewed Fedorov’s method 
and discussed critical angle phenomena for the elastic waves at an interface. These works 
mainly focussed on finding reflection and refraction angles using geometrical plotting on 
slowness surfaces. 
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More recently, RokWin et al [10] have described a unified approach to the 
numerical solution of the reflection-refraction problem for generally anisotropic media. ' 

Marker et al [11] presented a full numerical solution of the wave equation and 
showed that ray tracing method gives similar and consistent results. 

Kak [12] reviewed developments in Ultrasonic Computerized Tomography and 
has presented the reconstruction of velocity and attenuation coefficient images in medical 
applications. The effects of refi:action were ignored as only soft tissue-structures were 
studied. Digital ray-tracing algorithms have been developed by many workers to account 
for refraction in ultrasound tomography. They have used Snell’s law, Fermat’s principle 
and the Eikonal equation in geometrical optics for it. Kak and Andersen [13] have 
presented a review of digital ray tracing in two-dimensional refractive index field. The 
various aspects of numerical implementation like discretization, numerical interpolation 
schemes, step-size and the use of smoothing windows have been dealt comprehensively 
in the above-mentioned work. 

Bold et al [14] have compared the different integration methods with regard to 
the three algorithms reported by Kak and Andersen. They have found that fourth order 
methods give a good combination of accuracy along with simplicity and speed. 
Evidently, accuracy increases with the increase in the order of integration and with the 
decrease in step-size. 

Wang et al [15] have developed a technique for ray tracing in non-homogeneous 
media (both isotropic and anisotropic). Their approach is based upon Fermat’s principle 
and can be readily incorporated into tomographic reconstruction algorithms. They have 
also investigated the effect of ray-bending on the accuracy of acoustic tomograms. 

Denis et al [16] have shown that the reconstruction quality can be improved 
significantly for refractive index variations of upto 10% using curved-ray methods. They 
have presented both simulation and experimental results to verify the reduction of 
refraction artifacts by this method. 

Rathore et al [17] have carried out the tomographic reconstruction of an object 
containing two defect geometries at a fixed distance. In their work, the contrast i.e. the 
relative difference in material velocity or refiractive index between the defect and base 
material, is varied to study the reconstruction of two defects separately. In a later work. 
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Rathore et al [18] have varied the distance between two defect geometries so as to 

f 

establish the minimum gap at which they can be resolved. 

1.3 Present Work 

In the present work, the tomographic reconstruction of square, circular and elliptical 
defects in isotropic materials, using simulated time of flight data has been carried out. 
The quality of the reconstruction has been studied by varying the contrast, i.e. the relative 
difference in material velocity or refractive index between the defect and the base 
material. In the present work, several rays are considered from a single source and the 
one, which reaches nearest to the receiver, has been considered. The materials selected 
for the purpose of investigation are (1) Steel with Aluminium insert, (2) Steel with Nickel 
insert and (3) Steel with brass insert. Snell’s law has been used to accoimt for the 
bending of the ultrasonic wave inside the medium. Edge detection has been applied to 
find out the orientation of the interface for refraction in the presence of a curved artifact. 
The relaxation parameter has also been varied to study its influence on the convergence 
rate of the iterative algorithm. Plots indicating the paths taken up by the ultrasonic rays 
after several iterations and tables furnishing the values of the normalized RMS errors for 
different values of the relaxation parameter have also been provided in the work. 

1.4 Thesis organization 

The various chapters in the thesis deal with the following aspects. 

• Chapter 2 presents the problem formulation in some detail, the tomographic 
algorithms that have been used for the reconstruction, the other tomographic 
algorithms that are in vogue and finally a discussion on the computer implementation 
of these algorithms. 

• Chapter 3 deals with the effect of ray bending caused by the presence of refraction 
artifacts. A somewhat comprehensive discussion on edge detection used here to find 
out the orientation of the interface for refraction in the presence of a curved artifact 
has also been incorporated in this chapter. 

• Chapter 4 gives the reconstructed slowness images for the three discussed 
combinations of the parent material and the insert material. It also discusses the effect 
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of relaxation parameter on the convergence rate of the tomographic algorithm used in 

I 

this work. Tables supplying normalized RMS errors for different values of the 
relaxation parameter and plots showing typical ray-paths after several iterations have 
been provided along with this chapter. 

• Chapter 5 finally concludes the present work and gives suggestions for future work. 
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Chapter 2 


Reconstruction Methods and Problem Formulation 

2.1 Introduction 

The applications of computerized tomography (CT) in today’s engineering world are 
manifold. It has got immense applications in the field of medicine as well. In 
conventional projection radiography, the result is obtained in the form of a shadowgraph 
requiring a good interpretation. CT, however, uses tomographic reconstruction to 
combine information from multiple projections. Thus, the CT information is displayed 
quantitatively as a reconstructed slice plane of the object and it can provide accurate 
quantitative measures of material characteristics in a small volume of the component. 

Although computerized tomography was originally conceived for X-rays, in 
recent times, its horizon has considerably broadened with its successful implementation 
for other energy sources such as ultrasonics, lasers and magnetic resonance etc. Among 
these, ultrasonic tomography has special applications in the field of non-destructive 
evaluation. Ultrasonic tomography can be defined as a technique to determine values of a 
spatially varying acoustic parameter of interest across a slice through an object. 

2.2 Reconstruction methods 

Tomography can be classified into (a) transform (b) series expansion and (c) optimization 
methods. Transform methods generally require a large number of projections for a 
meaningful answer. However, in most of the cases it is not practicable to record such a 
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large number of projections partly owing to inconvenience and partly due to time and 
cost. Hence, as a rule requirement of a large number of projections is not desirable and 
one must look for methods that converge within a , few iterations. Limited-view 
tomography is best accomplished using the series-expansion method. As limited data 
tomography does not have a unique solution, the algorithms are expected to be sensitive 
to the initial guess of the field that start the iterations. Optimization-based algorithms are 
known to be independent of initial guess, but the choice of the optimization function 
plays an important role in the result obtained. Depending on the mathematical definition 
used, the entropy extremization route may yield good results, while the entropy 
minimization principle may be suitable in other applications. For the algebraic techniques 
considered in the present study, an unbiased initial guess such as a constant profile was 
seen to be good enough to predict the correct field variable. Tomography, being an 
inverse technique was seen to preserve the noise in the initial data. However, the 
dominant trend in the field variable was captured during tomographic inversion. Actually 
the problem of reconstruction is a problem of inversion of a rectangular matrix. Iterative 
techniques that are used in tomography can be viewed as developing a generalized 
inverse of the matrix. However, the matrices are mainly sparse ones. Thus general 
purpose matrix libraries can not be used to invert them sine they are highly ill- 
conditioned and rectangular in structure. The tomographic algorithms can be seen as a 
systematic route towards a meaningful inversion of the matrix. 

Series expansion methods have been adopted as the appropriate technique for 
image reconstruction. These methods are iterative in nature and consist necessarily of 
four major steps, namely: 

1 . Initial assumption of the field to be reconstructed 

2. Calculation of the correction 

3. Application of the correction 

4. Test for convergence 
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2.2.1 Algebraic reconstruction techniques 


source 



Figure 2.1, Intercept made by the i* ray in the j* pixel 


In Figure 2.1 a square grid on the image g(x,y) has been superimposed with the premise 
that in each element of the grid g(x,y) is constant.The field value in the m‘*’ cell is denoted 
by Xm . The total number of cells is taken to be N. In algebraic reconstruction techniques a 
ray is defined somewhat differently. In this context, a ray is now defined as a “fat” line 
running through the xy plane. The other important assumption behind the theory of 
algebraic reconstruction techniques is that the ray-width is approximately equal to the 
image cell width. In case of defect reconstruction for non-destructive evaluation 
purposes, the projections are basically the time-of-flight data i.e. the transit time between 
the source and the receiver. The relationship connecting the field value (xi) with the 
projection value (pO can be expressed as 

WyXj = p,,i = 1,2, ,M (2.1) 

y=i 

where “M” is the total number of rays (in all projections), and Wy is the weighting 
factor that represents the contribution of the j^ cell to the i* ray integral. Also “N” is the 
total number of cells in the grid. 
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Conventional matrix methods can be used to invert the system of equations in 
Eq.(2.1) for small values of “M” and “N”. However, in practice the values of “N” are 
quite large, sometimes of the order of 65,000 (for 256x256 images), and in most cases for 
images of the above-mentioned size, “M” also has the same magnitude. Thus the weight 
matrix becomes too large, which precludes any possibility of direct matrix inversion. 
Also the weight matrix is a sparse one with many of its elements being zero. General- 
purpose matrix libraries cannot be used to invert such matrices since they are highly ill- 
conditioned and rectangular in structure. Of course, when noise is present in the 
measurement data and when M<N, even for small N, it is not possible to use direct matrix 
inversion, and some least-squares methods may have to be used. But such methods are 
impractical for large values of both M and N. 

For large values of M and N, there exist very attractive iterative methods for 
solving the system of equations given in Eq. (2.1). The earliest methods were proposed 
by Kacmarz [19], which is commonly referred to, as the method of projections. Writing 
Eq. (2.1) in an expanded form, the computational steps involved in the method are 


explained. 

^13^3 ~ P\ ( 2 - 2 ) 

■*■^ 23^3 ~ P 2 (^- 3 ) 


'^^U2^2 '^^M3^3 ~ Pm (2-4) 


A grid representation with N cells gives an image N degrees of freedom. Thus, an image 
which is represented as (xi, X 2 ,.... xn) may be considered as a single point in an N- 
dimensional space. In this space each of the above equations represents a hyperplane. 
When a unique solution to this system of equations exists, the intersection of all these 
hyperplanes is a single point giving the solution. In this context for the purpose of 
explanation the case of two variables is considered. 

They satisfy the following equations: 



+ =Pi 

(2.5) 

1^21^1 

+ ^ 22^2 =P2 

(2.6) 
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The computational procedure for locating the solution to the system of equations consists 
of starting with an initial guess, projecting this initial guess onto the first line, 
reprojecting the resultant point on the second line and then projecting back onto the first 
line, and so on until the two successive solutions coincide. If a unique solution exists, the 
iterations will always converge to that point. 

This algebraic reconstruction algorithm generally called as algebraic reconstruction 
technique (ART), may further be classified into different categories depending upon the 
correction applied. Four classes of ART are mentioned below. 

1 . Mayinger ART 

2. Gordon ART. 

3. Gilbert ART 

4. Anderson ART 

In this problem, Mayinger ART has been used. In travelling from the source to the 
receiver, the ultrasonic ray can travel in a straight line or it can take a curved path 
depending on the difference in its velocities through the parent material and the insert 
material. If this difference exceeds 15%, taking recourse to ray bending gives better 
results than straight ray approach. In this problem ray-bending has been taken into 
consideration while tracing the path of the ultrasonic beam firom the source to the 
receiver. Chapter 3 deals with that phenomenon. 

MART 

There are various variants of the ART algorithms depending on the way the corrections 
are applied. One such modified version is the MART, which is also known as the 
multiplicative ART. In MART, after every iteration, a new estimate of the field value is 
calculated by multiplying the previous value with a correction term, i.e. 

7 = 1,2, 3, N (2.7) 

where “C” is the correction term, and is the slowness value of the j*'’ pixel at the end 
of the k* iteration. Again depending on the correction term, MART may be further 
classified into four classes. The underlying algorithm remains the same but the only 
difference lies in the correction term. The four algorithms along with the associated 
correction terms are provided below. 
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l.GBHMART 


Gordon, Bender and Herman proposed this algorithm for tomographic applications. In 
this algorithm, the correction term is given by 


C = l- 


2 


(K) 


ij / max 


1- 


isX 


^0,andW. ^0 


( 2 . 8 ) 


else 


C=1 


where, 

X is the relaxation parameter, 

th th 

Wij is the intercept made by the i ray in the j pixel. 

(Wij)max is the maximum of all intercepts made by the i'*’ ray. 

(xi)ni and (Ti)c are the measured and calculated time of flight data respectively. 


2. GH MART 

This algorithm was proposed by Gordon and Herman. The correction term is given by, 


C = l- 




A(f, 


1 - 




(2.9) 


else 


C=1 


3. LENT MART 


Lent has proposed another multiplicative updating scheme, in which the correction term 
is computed as. 


C = 


(lX 

(t), 






/ 


( 2 . 10 ) 


else 


C=1 

4. LENT MART2 

A slight modification of the above algorithm leads to Lent MART2. Here the correction 
term is given as. 
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( 2 . 11 ) 


else 





C=1 


where 

Wij* is given by Wi/Wmax 


SIRT 

The ART algorithms enjoy a rapid convergence in the root-mean-squared error criterion. 
However, the reconstructed images exhibit a very noisy salt-and-pepper characteristic. 
Smoother images are obtained by the simultaneous iterative reconstruction technique 
(SIRT) algorithm, in which the corrections in Eq.(2.16) are computed due to all the 
equations before altering any pixel values. Only then are the pixel values updated by their 
average corrections. This constitutes one iteration. In contrast with ART reconstructions, 
SIRT methods yield smoother images, although at the expense of slower convergence. 

In ART algorithms, the salt-and-pepper noise is caused by the inconsistencies 
introduced in the set of equations by the approximations commonly used for Wjk. This 
results in computed ray sums in Eq.(2.14) being a poor approximation to the 
corresponding line integrals. The effect caused by such inconsistencies is exacerbated by 
the fact that as each equation corresponding to a ray in a projection is taken up, it changes 
some of the pixels just altered by the preceding equation in the same projection. The 
SIRT algorithm also suffers from these inconsistencies in the forward process (appearing 
in the computation of qi in Eq.(2.14); but by the elimination of the continual pixel update 
as each new equation is taken up, smoother reconstructions can be obtained. 


2.3 Computer Implementation of the procedure 

Stepl 

An initial guess to the solution of the discussed system of equations is assumed. This 
guess, denoted by xi^°\ X 2 ^°\ ..., xn^°\ is represented vectorially by in the N- 
dimensional space. This initial guess is projected on the hyperplane represented by the 
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first of the equations i.e. Eq (2.2) thus yielding Then is projected on the 
hyperplane represented by the second equation to yield and this process is repeated 
till the entire set of equations is exhausted. When x^"’^ is projected on the hyperplane 
represented by the j* equation to yield x^\ the process can be mathematically described 
by 


= 


WjWj 


( 2 . 12 ) 


where Wj = (wji, wja wjn). 

However, this equation cannot be used in its present form. The difficulty with 
using the above equation lies with the problem it poses in the calculation, storage and fast 
retrieval of the weight coefficients wy. For reconstructing an image on a 128x128 grid 
from 150 projections with 150 rays in each projection, the number of elements in the 
weight matrix becomes enormously large which can pose major problems in fast storage 
and retrieval in applications in which reconstruction speed is important. 

To circumvent the implementation difficulties caused by the weight coefficients, 
recourse can be taken of different algebraic approaches. Many of these are 
approximations of Eq.(2.12).For a discussion of these approaches, Eq.(2.12) is first recast 
into the following equation. 





(2.13) 


where , 


( 2 . 14 ) 

k=l 

The underlying logic behind these equations is that, on projecting the (j-l)* solution onto 
the j* hyperplane, the pixel value of the m* element in the grid, whose current value is 
, is obtained by correcting this current value by A Xm^^ , where 
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(2.15) 


k=} 

Here, as defined earlier, pj is the measured ray-sum along the j* ray. Similarly, qj can 
be considered to be the computed ray-sum for the same ray based on the (j-1)* solution 
for the image pixel values. The correction A Xm to the m* cell is obtained by first 
calculating the difference between the measured ray-sum and the computed ray-sum, 
normalizing this difference by Ewjk^, and then assigning this value to all the image cells 
in the ray, each assignment being weighted by the corresponding w*"’. 

In one approximation to Eq.(2.15), the wjk are simply replaced by I’s and O’s, 
depending on whether the center of the k'*' image cell is within the j* ray. This makes the 
implementation easier, because such a decision can easily be made at computer run-time. 
Clearly in this case, the denominator in Eq.(2.15) can be replaced by Nj, which is the 
number of image cells whose centers are within the ray. The correction to the m* 
image cell from the j equation in Eq.(2. 1 5) can now be written as 


Ax^ = 


(P i -^/) 




(2.16) 


for all the cells whose centers are within the j*** ray. The approximation in Eq.(2.16), 
although easy to implement, may lead to artifacts in the reconstructed images if Nj is not 
a good approximation to the denominator of Eq.(2.15). 


Step 2 

The process of taking projections on different hyperplanes is continued until is 
obtained, which in turn is obtained by taking the projection on the last of the system of 
equations i.e. Eq.(2.4). This completes one cycle of this iterative procedure. 

In steps 1 and 2, as each new equation starting from Eq.(2.2) and ending at 
Eq.(2.4) is taken up, the values of all the affected pixels are updated.In a related 
approach, which may lead to fewer artifacts, the changes in Axm^^ caused by the j**’ 
equation in the system of equations, can again be computed. However, now the value of 
the m* cell is not changed at this time. The cell values remain intact until all the 
equations are finished and then only at the end of each iteration the cell values are 
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changed, the change for each cell being the average value of all the computed changes for 
that cell. ' 

Step 3 

Iterations are continued by projecting on the first hyperplane again. This process 
continues till all the hyperplanes have again been cycled through, resulting in The 
second iteration is started by projecting on to the first hyperplane, and so on. 
Iterations are stopped when the computed changes of the image cell values are negligible 
firactions of their current values. 

2.4 Formulation of the present problem in the light of the algebraic 
reconstruction algorithm 

The experiment consists of generating ultrasound and thereby sending the same through 
an isotropic material. The positions of the source and receiver are first fixed up according 
to the angle, which the ultrasonic ray takes in travelling through the medium. As already 
mentioned earlier in this chapter, a fan-beam configuration has been taken with five rays 
emanating from each point-source. Finally among these five rays, only one ray, which 
reaches nearest to the corresponding receiver position, is considered and the other four 
rays are discarded. Thus in the final analysis, one ray has been considered for a source- 
receiver pair. The algebraic reconstruction algorithm, which has been discussed in the 
previous sections of this chapter, is applied to reconstract the slowness field in a slice- 
plane of the isotropic material. For isotropic materials, slowness is simply the inverse of 
the longitudinal sonic velocity through the medium. The projection data for this problem 
is the time of flight, i.e. the transit time between the source and the receiver. The iterative 
process starts from an initial approximation of slowness, S° , and the current iterate is 
corrected to a new iterate by taking into account only a single ray and changing the 
image values of the pixels through which the ray passes. The disparity between the 
measured time of flight, and the pseudo-projection data, also called the numerical time of 
flight, obtained from the current image S'^, is redistributed among the pixels falling in the 
path of the ray, proportionally to their weights, which are the intercepts that the ray 
makes with the pixels through which it passes. In this way, the pixel values along a 
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particular ray are corrected to conform to that measurement of the projection data without 
changing the rest of the image. The whole process is described mathematically as: 

z-< > 


+■ 


w 


‘ l|2 


■W 


(2.17) 


where, 

w^±w, 

and 


(2.18) 

(2.19) 


II w' f > (2.20) 

Introducing the relaxation parameter, the modified step for updating the iterate can be 
expressed as 

( 2 . 21 ) 


where, % is the relaxation parameter, which is a real number, usually lying between 0 and 
- 2 . 

If the correction is negative, it may happen that the calculated field of a 
reconstructed image is negative, in which case it is set to zero, as discussed by Gordon. 


2.5 Closure 

A description of the commonly applied tomographic algorithms has been presented in 
this chapter. The problem at hand has also been formulated in the light of these 
algorithms. The various steps involved in the computer implementation of these 
algorithms are also provided in this chapter. 
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Chapter 3 


Ray-Bending and Edge-Detection 

3.1 Introduction 

The ultrasonic wave undergoes refraction when it encounters an inhomogeneity in its 
path of propagation. This phenomenon becomes more complex if there are more than one 
inhomogeneity present in a solid object. The use of digital ray tracing method has paved 
the way for the determination of the actual path of wave propagation between the source 
and the receiver under such condition. Using digital ray tracing method in conjunction 
with the algebraic reconstruction algorithms, it is possible to reconstruct either the 
slowness or the velocity field in a material with refraction artifacts. 

3.2 Ultrasonic wave propagation 

Ultrasonic waves of small amplitude propagate as linear elastic waves whose propagation 
through a medium results in particle oscillations about their equilibrium positions. In 
ideal liquids and gases only one mode of propagation exists. This longitudinal mode has 
particle displacements parallel to the direction of ultrasonic propagation. The ultrasonic 
velocity is governed by the compressibility of the materials and for these materials, is 
independent of propagation direction. Elastic solids can accommodate transverse waves 
(where the particle displacements are perpendicular to the direction of propagation), as 
well as longitudinal waves. The particle motions in elastic solids associated with these 
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bulk inodes of propagation can be resolved into three perpendicular components, one 

I 

longitudinal and two transverse to the direction of propagation. 

Wave propagation in anisotropic solids 

In anisotropic elastic solids, each of the three modes, one longitudinal and two transverse 
modes has its own characteristic velocity, although in isotropic solids the two transverse 
velocities degenerate into one. For anisotropic solids, the phase velocity (v) of the three 
possible bulk modes (quasi-longitudinal, fast and slow quasi-shear) in a direction defined 
by the unit propagation vector, n, with components nj, j = 1,2,3 can be found from the 
eigenvalues of the Christoffel equation given as follows: 

)ij = Oi = 1.2,3. (3.1) 

where p is the density of the material, k = w/v, is the wave-number, w is the angular 
frequency, cyki are the components of the elastic stiffness tensor, dj, j = 1,2,3 are the 
components of the unit particle displacement vector, d, and by is the Kronecker delta. For 
a given propagation direction vector n, the squares of the phase velocities of the bulk 
modes are obtained as the eigenvalues of Eq.(3.1). The particle displacement vector, d, of 
a bulk mode is obtained as a unit eigenvector for the phase velocity (eigenvalue) of the 
bulk mode. The slowness (or inverse velocity) surface gives the inverse phase velocity as 
a function of propagation direction and is independent of the angular frequency. 

An important consequence of anisotropic wave propagation is that the energy of 
a wave packet does not necessarily propagate parallel to the direction of wave 
propagation, unlike in an isotropic media, where both directions are coincident. In a 
lossless medium, the energy or group velocity is perpendicular tp the slowness surface. 
The powerflow angle defined as the angle between the energy velocity and the 
propagation direction depends on the shape of the slowness surface. The group velocity, 
which governs the speed of a modulated disturbance (i.e. a pulse), can be determined 
once the phase velocity and the displacement vector for a specific mode are known. 
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3.3 Reflection and transmission at an interface 


When an ultrasonic wave propagating in a homogeneous medium passes through an 
interface at a normal incident angle where the velocity (i.e. the elastic constant and 
density) of the second medium is different (i.e. solid-solid or solid-liquid interface), the 
incident wave is reflected and transmitted without any change in the propagation 
direction. However, when the ultrasonic wave, propagating in a homogeneous medium 
passes through an interface at a non-normal incident angle, the incident wave is mode 
converted (a change in the nature of the wave motion) and refracted (a change in the 
direction of wave propagation). These phenomena may affect the entire beam or only a 
small portion of the beam and the sum total of the changes that occur at the interface 
depends on the angle of incidence and the velocity of the ultrasonic waves leaving the 
point of impingement on the interface. The reflection and refraction of these ultrasonic 
waves is governed by the acoustic analog of Snell’s law (the ratio of the sine of the angle 
of incidence to the sine of the angle of reflection or refraction equals the ratio of the 
corresponding wave velocities). If two media are represented by 1 and 2 , then Snell’s law 
can be put in mathematical form as follows: 


sin 6*1 _ (sin^;), _ (sin 6*2)^ 

(^2)/ (^2), 


where 61, (62)1 and (02)s are the ray angles with respect to the surface normal and vi_ (V2)i 
and (V2)s are the velocities of the longitudinal wave in the first medium, the longitudinal 
and shear wave velocities in the second medium. The phenomenon has been pictorially 
depicted in Fig.( 3 . 1 ). However, Snell’s law is equally valid for a solid-solid interface. For 
the problem at hand, a solid-solid interface is present. If the angle of incidence is small, 
the incident wave undergoes a mode conversion at the interface, resulting in the 
simultaneous propagation of longitudinal and shear waves in the second medium. If the 
angle of incidence is increased, the direction of the refracted longitudinal wave will 
approach the plane of the interface. At the first critical incidence angle, the refracted 
longitudinal wave will disappear leaving only a refracted (mode converted) shear wave to 
propagate in the second medium. Now if the incidence angle is increased beyond this 
critical value, the direction of this refracted shear wave will approach the plane of the 
interface. At this second critical incidence angle, the refracted shear wave will disappear 
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and the incident wave is only reflected at an angle equal to the incident angle. This 

I 

methodology applies to any interface, i.e. liquid-solid, solid-solid and solid-liquid. 



Figure (3.1) Reflected and refracted rays at a solid-solid interface 


3.4 Ultrasonic ray theory 

Geometrical acoustics describe the propagation behaviour of acoustic rays in liquid and 
solid bodies. In the development of ultrasonic sensor methodology for reliable image 
reconstruction using refracted ray paths, the sometimes severe nature of ray bending 
implies a need for precise knowledge of ray paths, particularly when time of flight data is 
needed. 

In general, there are two methods to determine the ray paths between two points; 
shooting and bending. The shooting method uses an iterative procedure to determine the 
ray path from a source point by solving the differential equations, that follow from ray 
theory for different initial conditions until the ray arrives at the receiver point. Ray 
bending uses Fermat’s principle as a starting point. Ray bending is assumed to occur at a 
refraction artifact. The ‘effect of ray bending cannot be ignored when the sonic velocity 
difference between the base and the insert material is over 15%. 

Both methods have certain limitations. By shooting a fan of rays from the 
source, an impression of the wave field can be obtained. However, convergence problems 
are known to occur, frequently in three dimensions. Also shooting will not find diffracted 
ray paths or ray paths in shadow zones, where ray theory breaks dovwi. With bending, 
every ray path satisfying Snell’s law and Fermat’s principle can be found, even a 
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diffracted one, but only for one source-receiver pair at a time. In the present work the 
refracted ray paths have been determined using Snell’ s law. 

3.5 Image Analysis and Edge Detection 

The reconstructed image obtained is refined using various edge detection techniques. 
Changes or discontinuities in an image amplitude attribute such as luminance or 
tristimulus value are fundamentally important primitive characteristics of an image 
because they often provide an indication of the physical extent of objects within the 
image. Local discontinuities in image luminance from one level to another are called 
luminance edges. In this section, several techniques for detecting the three basic t 3 qDes of 
discontinuities in a digital image: (1) points (2) lines and (3) edges will be discussed. The 
most common way to look for discontinuities is to run a mask through the image. This 
procedure involves computing the sum of products of the coefficients with the pixel 
values contained in the region encompassed by the mask. Thus, the response of the mask 
at any point in the image is given as follows: 

R = WiZ, + WjZj + ■'^9^9 (3-3) 

where Zi is the value of the pixel associated with the 3x3 mask coefficient Wj . The 
response of the mask is defined with respect to its center location. When the mask is 
centered on a boundary pixel, the response is computed by using the appropriate partial 
neighbourhood 

3.5.1 Point detection 

The detection of isolated points in an image is straightforward. 'A point is said to be 
detected at the location on which a mask is centered if |R1 > T, where T is a nonnegative 
threshold, and R is given by Eq.(3.3). Thus the imderlying idea here is to compute the 
weighted differences between the center point and its neighbours. The idea is that the 
pixel value of an isolated point will be quite different from the pixel value of its 
neighbours. 
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3.5.2 Line detection 


The next level of complexity involves the detection of lines in an image. The 
phenomenon is best demonstrated by considering the four simple masks given below. 


-1 

-1 

-1 

2 

2 

2 

-1 

-1 

-1 


Horizontal 


-1 

-1 

2 

-1 

2 

-1 

2 

-1 

-1 


+45 degree 



If the first mask were moved around an image, it would respond more strongly to lines 
oriented horizontally. With constant background, the maximum response would result 
when the line passed through the middle row of the mask. A similar logic would reveal 
that the second of the masks shown above responds best to lines oriented at 45°; the third 
mask to vertical lines; and the fourth mask to lines in the -45° direction. These directions 
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3.5.2 Line detection 


The next level of complexity involves the detection of lines in an image. The 
phenomenon is best demonstrated by considering the four simple masks given below. 




If the first mask were moved aroimd an image, it would respond more strongly to lines 
oriented horizontally. With constant background, the maximum response would result 
when the line passed through the middle row of the mask. A similar logic would reveal 
that the second of the masks shown above responds best to lines oriented at 45°; the third 
mask to vertical lines; and the fourth mask to lines in the -45° direction. These directions 
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can be easily established by noting that the preferred direction of each mask is weighted 
with a larger coefficient (i.e. 2) than other possible directions. 

Let Ri, R 2 , R 3 and R 4 denote the responses of the masks shown above, where R 
has been defined in Eq.(3.3). For example, after running all the masks through an image 
if at a certain point in the image, | Ri | > | Rj |, for j = 2,3,4 that particular point is said to 
be more likely associated with a horizontal line. 

3.5.3 Edge detection 

Although point and line detection certainly are important elements of image 
segmentation, edge detection is by far the most common approach for detecting 
meaningful discontinuities in an image represented by pixel values (gray levels). The 
reason is that isolated points and thin lines are not frequent occurrences in most practical 
applications. 

An edge is the boundary between two regions with relatively distinct gray-level 
properties. The commonly applied segmentation techniques are valid on the assumption 
that the regions in question are sufficiently homogeneous so that the transition between 
two regions can be determined on the basis of gray-level discontinuities alone. 

Basically, the idea underlying most edge-detection techniques is the computation 
of a local derivative operator. Actually, the magnitude of the first derivative can be used 
to detect the presence of an edge in the image and the sign of the second derivative can be 
used to determine whether an edge pixel lies on the dark or light side of an edge. The 
second derivative has always a zero crossing at the midpoint of a transition in gray level. 
The zero crossings provide a powerful approach for locating edges in an image. This is in 
context to a one-dimensional horizontal profile. However, a similar argument applies to 
an edge of any orientation in an image. A profile perpendicular to the edge direction at 
any desired point can be defined always and the result can be interpreted in similar lines 
to the horizontal unidimensional profile. The first derivative at any point in an image is 
obtained by using the magnitude of the gradient at that point. The second derivative is 
similarly obtained by using the Laplacian. 
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First order derivative edge detection 

There are two fundamental methods for generating first order derivative edge gradients. 
One method involves generation of gradients in two orthogonal directions in an image, 
while the other utilizes a set of directional derivatives. 

Orthogonal gradient generation 

An edge in a continuous domain edge segment F(x,y), can be detected by forming the 
continuous one-dimensional gradient G(Xjy) along a line normal to the edge slope, which 
is at an angle 0 with respect to the horizontal axis. If the gradient is sufficiently large, i.e. 
above some threshold value, an edge is deemed present. The gradient along the line 
normal to the edge slope can be computed in terms of the derivatives along orthogonal 
axes according to the following equation as follows: 

= (3.4) 

Ax Ay 

Thus an edge gradient G(j,k) can be written in terms of a row edge gradient GR(j,k) and a 
colunrn gradient GcO ,k). The spatial gradient amplitude is given by 

|[G,0'X)]^ + [Gc0'.i)f ( 3 - 5 ) 


For computational efficiency the gradient amplitude is sometimes approximated by the 
magnitude combination 

G(j, k) = (y, k} + |Gc (;, k} (3.6) 

The orientation of the spatial gradient with respect to the row axis is 

The remaining issue for discrete domain orthogonal gradient generation is to choose an 
accurate discrete approximation to the continuous differentials along the horizontal and 
vertical directions. 
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The simplest method of discrete gradient generation is to form the running 
difference of pixels along rows and columns of the image. The row gradient and column 
gradient are defined as 


GR{j\k) = F{j,k)~F{j,k-\) (3.8) 

Gc{j,k)=F(j,k)-F{j + \,k) (3.9) 

These definitions of row and column gradients, and subsequent extensions are so chosen 
such that Gr and Gc are positive for an edge that increases in amplitude from lefl-to-right 
and from bottom-to-top in an image. 

Diagonal edge gradients can be obtained by forming running differences of 
diagonal pairs of pixels. This is the basis of the Roberts cross difference operator, which 
is defined in magnitude form as 

G{j,k) = \G,{j,k\ + \G^(j,k} (3.10) 

and in square root form as 

where 

G,{j,k)=F{j,k)-F{j + \,k + \) (3.12) 

Gj (;, k) = F(j, k + \)-F{j + U) (3.13) 


The edge orientation with respect to the row axis is 


e{j,k) = ^ + tan' 




(3.14) 


Visually, the objects in an image appear to be better distinguished with the Roberts 
square root gradient than with the magnitude gradient. 

Prewitt has introduced a 3x3 pixel edge gradient operator. The Prewitt square 
root edge gradient is defined as 
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with 


G(j,k)^\^G,(j\k)f + [a,lj,k)fV'^ (3.15) 

G, [j, [(A + H + A )- (A + ^7 + A)] (3-16) 

^cOA) = j^^[(A +^i + A)~(A +^5 + a)] (3-17) 

where k = 1 and Ai are given as follows. 


Ao 

Ai 

A 2 

A7 

FO^k) 

A 3 

A 6 

A, 

A 4 


Numbering convention for 3x3 edge detection operators. 

In this formulation the row and column gradients are normalized to provide unit gain 
positive weighted and unit gain negative weighted averages about a separated edge 
position. 

The Sobel edge detector, which is used in the present problem, differs from the 
Prewitt edge detector in that the values of the north, south, east and west pixels are 
doubled, i.e. k=2. The motivation for this weighting is to give equal importance to each 
pixel in terms of its contribution to the spatial gradient. Another idea is to give north, 
south, east and west weightings a value of 1.414 so that the gradient is the same for 
horizontal, vertical and diagonal edges. The Prewitt operator is more sensitive to 
horizontal and vertical edges than diagonal edges; the reverse is true for the Sobel 
operator. 

A limitation common to the edge detection operators is their inability to detect 
accurately edges in high-noise environments. This problem can be alleviated by properly 
extending the size of the neighbourhoods over which the differential gradients are 
computed. 
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3.5.4 Threshold selection 

After the edge gradient is performed for the differential edge detection methods, the 
gradient is compared to a threshold to determine if an edge exists. The threshold value 
determines the sensitivity of the edge detector. For noise-free images the threshold can be 
chosen such that all amplitude discontinuities of a minimum contrast level are detected as 
edges and all others are called non-edges. With noisy images, threshold selection 
becomes a tradeoff between missing valid edges and designating noise-induced false 
edges. 

3.5.5 Second order derivative edge detection 

Second order derivative edge detection techniques employ some form of spatial second 
order differentiation to accentuate edges. An edge is marked if a significant spatial 
change occurs in the second derivative. Two types of second order derivative methods are 
generally used: (1) Laplacian and (2) Directed second derivative. 


Laplacian generation 

The edge Laplacian of an image function F(x,y) in the continuous domain is defined as 


■S^F{x,y) -5^F{x,y) 




(3.18) 


The Laplacian G(x,y) is zero if F(x,y) is constant or changing linearly in amplitude. If the 
rate of change of F(x,y) is greater than linear, G(x,y) exhibits a sign change at the point of 
inflexion of F(x,y). The zero crossing of G(x,y) indicates the presence of an edge. 

In the discrete domain the simplest approximation to the continuous Laplacian is 
to compute the difference of slopes along each axis. 

G{j,k) = [fO’, A:)- F{j,k-l)]-[F(j,k + 1)- /^0-,A:)]+ + 1,^)]- [^0' 

(3.19) 

In the present problem, the above approximation has been taken for the continuous 
Laplacian. 
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Directed second order derivative generation 

Laplacian edge detection techniques employ rotationally invariant second order 
differentiation to determine the existence of an edge. The direction of the edge can be 
ascertained during the zero crossing detection process. An alternate approach is to first 
estimate the edge direction, and then compute the one-dimensional second order 
derivative along the edge direction. A zero crossing of the second order derivative 
specifies an edge. 

The directed second order derivative of a continuous domain image F(x,y) along 
a line at an angle a with respect to the horizontal axis is given by 




Unlike the Laplacian, the directed second order derivative is a nonlinear operator. A key 
factor in the utilization of the directed second order derivative edge detection method is 
the ability to accurately determine its suspected edge direction. 


3.6 Closure 

In this chapter ray bending has been discussed. Snell’s law, which dictates the direction 
taken up by the ray when it encounters an inhomogeneity in its path of propagation, has 
also been taken up. Finally various image segmentation techniques, namely, point 
detection, line detection and most importantly edge detection has been discussed at 
length. The edge detection techniques, which have been used in this present work, have 
also been mentioned. 
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Chapter 4 


Results and discussion 

4.1 Introduction 

In the present work, Mayinger ART algorithm has been used for the reconstruction of the 
tomographic images. Here essentially the slowness (inverse of ultrasonic velocity) 
distribution in the domain from numerous projection data are computed. The projection 
data basically consists of measuring the time of flight from the source to the receiver. 
Curved ray approach has been used to compute the ray paths. The acoustic analogue of 
Snell’s law has been used to ascertain the path taken by the ray after it encounters a 
refraction artifact. 

The present study has been conducted for three different combinations of parent 
material and insert materials. They are (1) Steel with Aluminium insert (2) Steel with 
Nickel insert and (3) Steel with brass insert. Three different t 3 ^es of defect geometries 
have been considered, namely, square, circular and elliptical. The quality of 
reconstruction has been studied with respect to the contrast of the base and insert 
materials, i.e. the relative difference in wave velocities between the defect and the base 
material. The quality of reconstruction has also been found for values of relaxation 
parameters lying within the range 0.1 to 0,4. All the reconstructed images have been 
shown for these values of the relaxation parameter only. The following table lists the 
properties of the materials, which have been used in this work. 
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Serial 

Number 

Base 

material 

Insert 

material 

Velocity of 
sound in 
base 
material 

Velocity of 
sound in 
insert 
material 

Projection 
data typp 

1 


Aluminium 

5160 m/sec 



2 

Steel 

Nickel 

5160 m/sec 

4820 m/sec 

4- view 

3 

Steel 

Brass 

5160 m/sec 

3700 m/sec 

4-view 


Table 4.1 Properties of the base material and the insert materials. 


4.2 Details of Test cases 


In the present work, the tomographic reconstruction of square, circular and elliptical 
defects have been carried out for three different combinations of parent material and the 
insert material. All the reconstructions have been carried out taking a 4- view projection. 











n 











■ 










■ 










■ 




























1 

i 










1 

1 

1 


1 

2 

3 






1 



Figure 4. 1 Grid for discretizing the domain 


The iterative reconstruction methods require the discretization of the domain to be 
represented by a square or rectangular grid. For the problem at hand, a 20 x 20 square 
grid hqs been taken. However, for clarity, a 10x10 grid has been shown in Figure 4.1. The 
individual elements in the grid are referred to as pixels. The pixels were numbered from 1 
to 400. Here the analysis is carried out by taking a four-view projection. The rays of 
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sound energy were sent through the sample at four different angles, namely, 0°, 45°, 90° 
and 135° with respect to the positive X-axis. While ultrasound traverses through the 
material, the rays pass through the various pixels. When the i* ray passes through the j**' 
pixel, the length of the intercept made by the i*^ ray in the j* pixel is denoted by wy, 
which is popularly called the weight. This has been depicted in Figure 4.2. In the present 
work, 20 rays each inclined at 0° and 90° have been taken along with 40 rays each 
inclined at 45° and 135° respectively. Thus in all, 120 rays have been considered for the 
image reconstruction. However, while collecting simulated time of flight data for four 
transmitter-receiver locations at 45° and 135°, the rays starting from those particular 
transmitter locations have not reached the corresponding receiver locations due to large 
refraction at the intervening defects. Thus, 1 12 rays have been considered finally. 



Figure 4.2. Intercept made by the i^ ray in the pixel 


In this problem, to trace the correct ray, a set of five rays emanating from each point 
source, with an angular spacing of 1.25° between two consecutive rays has been taken. 
Here are the details of the rays along each view. 
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1, Ray-Paths along 0° direction (x-axis) 

In this view, in all 20 data have been collected from 20 ray-paths, each source being 
located inside a pixel situated in the first column of the grid. The spacing between the 
rays of two consecutive sources has been taken as unity. The source points act as the 
transmitter points, while the corresponding receiver points are situated exactly opposite to 
the transmitter points. A set of five individual rays are sent out from a source at a spatial 
separation of 1.25°. Typical ray-paths for three 0° fan-beams, one passing through the 
Nickel insert, the other two just lying on either side of the insert are shown in Figure 4.4. 
The insert is circular in shape and the plot of rays after 5 iterations is shown in the figure. 
From the plot it can be observed that the central fan-beam passes through the insert where 
there is a marked deviation of the rays. Inside the insert, the rays are more closely spaced 
than in the region outside. 

2. Ray-Paths along 90° direction (y-axis) 

In this view also, in all 20 data have been collected from 20 ray-paths, each source being 
located inside a pixel situated in the first row of the grid. Here also the spacing between 
the rays of two consecutive sources has been taken to be unity. The source points act as 
the transmitter points, while the corresponding receiver points are situated exactly 
opposite to the transmitter points. A set of five individual rays are sent out from a source 
at a spatial separation of 1.25°. Typical ray-paths for three 90° fan-beams, one passing 
through the Aluminium insert, the other two just lying on either side of the insert are 
shown in Figure 4.3. The insert is circular in shape and the plot of rays after 5 iterations is 
shown in the figure. Here inside the insert, the rays are not so closely spaced as in the 
previous case. 


3. Ray-Paths along 45° direction 

For this configuration 40 source points have been considered. Here the transmitter points 
or the source points as they are referred to here, are located both on the left vertical and 
bottom horizontal edges of the scanned specimen. The corresponding receiver locations 
are on the top horizontal and right vertical edges of the sample. The spacing between the 
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rays of two consecutive sources is now 0.7071 units. A set of five individual rays are sent 

f 

out from a source at a spatial separation of 1.25°. Out of these five rays, the ray which 
reaches closest to the corresponding receiver position is taken and the rest are discarded. 
However, while collecting simulated time of flight data for four transmitter-receiver 
locations, the rays starting from those particular transmitter locations have not reached 
the corresponding receiver locations due to large refraction at the intervening defects. 
Thus, 36 rays have been considered finally. Typical ray-paths for three 45° fan-beams, 
one passing through the Brass insert, the other two just lying on either side of the insert 
are shown in Figure 4.5. From the figure, it can be clearly observed that the central fan- 
beam passes through the insert for which there is marked deviation for all the rays. The 
insert is circular in shape and the plot of rays after 5 iterations is shown in the figure. 


4. Ray-paths along 135° direction 

For this configuration also, 40 source points have been considered. Here the transmitter 
points or the source points as they are referred to here, are located both on the right 
vertical and bottom horizontal edges of the scanned specimen. The corresponding 
receiver locations are on the top horizontal and left vertical edges of the sample. Here 
also, the spacing between the rays of two consecutive sources is 0.7071 units. A set of 
five individual rays are sent out from a source at a spatial separation of 1.25°. Out of 
these five rays, the one which reaches closest to the corresponding receiver position is 
taken and the rest are discarded. Here also, due to the absence of time-of-flight data for 
four positions of the source point, in the final analysis only 36 data have been taken into 
consideration. 

4.3 Image Reconstruction 

In Figures 4.3-4. 5, typical ray-paths for three sets of rays, one at 0°, the other at 90° and 
the third one at 45° are shown. One set is so chosen that it passes through the insert and 
the other two lie just on either side of the insert. The insert is circular in shape. There is a 
clear deviation of the rays when they pass through the insert indicating its presence in the 
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base material. The deviations are marked for the fan-beam inclined at 45°. All the plots 
have been taken after 5 iterations. Besides that, in Figures 4.6 and 4.7, the entire ray- 
paths are plotted, one for circular Aluminium insert and the other for circular Nickel 
insert. The plots are taken after 10 iterations. From Figures 4.3-4.5, it can be clearly 
observed that the central set of rays is so chosen that it entirely passes through the 
circular insert. Also in all the three cases, the deviation suffered by all the rays 
constituting the central set is the maximum, thus clearly indicating the presence of the 
insert in that region. Some of the typical slowness values obteiined after a few iterations 
are shown in section 4.4. Though the order of the slowness matrix is 20x20, only a 10x10 
submatrix has been shown here. 

4.3.1 Images obtained without ray bending consideration 

Several images obtained without using curved ray approach have been shown in Figures 
4.10-4.12 for the purpose of comparison with images obtained using ray bending. The 
images have been shown for all the three defect geometries taken. The difference in the 
quality of reconstruction is about 3% in terms of the normalized RMS errors for the 
Aluminium and Nickel inserts. In other words, this difference is quite acceptable because 
of the small difference in longitudinal sonic velocities between them and the base 
material i.e. Steel. However, the difference is much higher, about 5% for the case of 
Brass insert because there is more than 20% velocity difference of sound through Steel 
and Brass. Thus for Brass, the reconstructed images are of inferior quality for all the 
types of defect geometries. 

4.3.2 Images obtained with ray bending consideration 

The images obtained using ray bending are much superior than the same obtained without 
using ray bending. This is proved by the quality of reconstructed images. Also the 
normalized RMS errors are consistently lower by about 3% for all the defect geometries 
considered. However, the errors are about 5% lower for the case of Brass insert. For 
Brass insert though, ray bending consideration fails to produce superior images because 
of the high difference in longitudinal sonic velocity (20%) between Steel and Brass. But 
the quality of reconstruction is definitely better with ray bending consideration for all the 
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insert geometries. All the reconstructed images have been shown for four values of the 
relaxation parameter, namely, 0.1, 0.2, 0.3 and 0.4. The reconstruction is found to be best 
for values of relaxation parameter lying within this range. The reconstructed images 
obtained are also compared with the contour plots obtained for all the three defect 
geometries. The contour plots agree quite satisfactorily with the reconstructed images. 
Figures 4.13-4.24 show the images obtained using Aluminium insert, the first four 
showing the images for the square insert for four values of the relaxation parameter, the 
next four showing the same for the circular insert and the last four for the elliptical insert. 
The contour plots have been immediately provided after that in Figures 4.25-4.27 for the 
purpose of comparison with the actual images. Figures 4.28-4.39 show the images for the 
Nickel insert, the first four for the square insert, the next four for the circular insert and 
ending with the elliptical insert. Figures 4.40-4.51 show the same for the Brass insert in 
the same order. 


4.4 Slowness values 


Some representative values of slowness are shown for the case of circular Aluminium 
insert in Steel. Though the order of the slowness matrix is 20x20, a 10x10 submatrix is 
shown here. These values are as follows: 


220.0725 242.2074 
231.8917 181.8477 
215.2587 245.4265 
225.0943 258.0464 
314.3125 275.8070 
301.7101 281.7370 
257.1631 319.7945 
203.7310 307.8354 
308.3214 245.9457 
282.0224 234.3150 


171.4287 156.5328 
217.9228 204.4597 
248.1114 251.5816 
275.3279 231.8135 
200.1192 165.4211 
212.8128 161.6447 
234.7083 185.6201 
264.7107 310.0605 
236.8223 274.0878 
199.1402 222.6359 


226.5117 247.9110 
222.1963 204.0383 
177.6984 170.6863 
193.9398 177.7762 
170.3677 178.5341 
176.8761 171.0957 
174.7852 167.1565 
193.8991 203.8077 
254.9149 244.1753 
283.1949 243.8067 


191.4956 153.7925 
204.4881 201.2777 
212.2036 234.5763 
194.0678 237.1187 
150.6968 134,6873 
140.3271 126.5760 
147.0800 134.9354 
238.4700 234.6395 
235.1453 238.4317 
160.8466 148.9963 


238.3172 212.6100 
188.7148 237.0013 
240.6320 228.3982 
207.2527 185.5086 
260.5496 233.4836 
151.0003 240.6375 
144.4117 184.5325 
164.0039 162.8761 
237.9478 231.5922 
232.9315 256.5932 


The values shown are obtained after 10 iterations. 
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4.5 RMS error 


The following tables give the values of the normalized RMS errors for different values of 
the relaxation parameter and for different combinations of the insert and the base 
material. The RMS error measures the deviation of the reconstructed field from the exact 
field. It is defined as: 

' 71 ) -(s ) 

max /a min Ja 

where, Sia and Sir are the actual and reconstructed slowness values of the pixel 
respectively. N is the total number of pixels. (Smax)a and (Smin)a are the actual maximum 
and minimum of all pixel values and Erms is the RMS error. 


Serial Number 

Relaxation 
Parameter (X) 

Number of 
iterations 

Normalized RMS 

error 

1 

0.100 

450 

18.136 

2 


425 

18.812 

3 



18.818 

4 

0.200 


18.761 

5 

0.250 

363 

18.323 

6 

0.280 

341 

18.519 

7 

0.300 

323 

18.131 

8 

0.345 


18.751 

9 


283 

18.121 


0.350 

250 


11 

0.371 

210 

18.161 

12 

0.380 

203 

18.321 

13 

0.400 

178 

18.112 

14 

0.500 

253 

18.891 

15 

1.000 

358 

19.121 

16 

1.500 

316 

19.361 

17 

1.800 

173 

19.034 


Table 4.2 showing normalized .RMS errors of the reconstructed images for square 
Aluminium insert implanted in steel from 4-view projection data with different relaxation 
parameters. 
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Serial Number 

Relaxation 
Parameter (X) 

Number of 
iterations 

Normalized RMS 

error 

1 

0.100 

437 

17.675 

2 

0.150 

417 

17.234 

3 

0.200 

356 

17.060 

4 

0.250 

323 

16.580 

5 

0.300 

265 

18.123 

6 

0.350 

234 

17.451 

7 

0.400 

179 

17.231 

8 

0.450 

145 


9 


132 


10 

0.550 

118 

17.941 

11 

0.600 

102 

17.894 

12 

0.650 

96 


13 

0.700 

89 

18.145 


Table 4.3 showing normalized RMS errors of the reconstructed images for circular 
Aluminium insert implanted in steel from 4-view projection data with different relaxation 
parameters. 


Serial Number 

Relaxation 
Parameter (X) 

Number of 
iterations 

Normalized RMS 

error 

1 

0.100 

421 

17.565 

2 

0.150 

411 

17.452 

3 

0.200 

373 

17.160 

4 

0.250 

306 

16.280 

5 

0.300 

279 

18.156 

6 

0.350 

251 

16.347 

7 

0.400 

192 

16.938 

8 

0.450 

163 

17.003 

9 

0.500 

146 

17.812 

10 

0.550 

131 

17.976 

11 

0.600 

112 

17.954 

12 

0.650 

104 

17.937 

13 

0.700 

98 

18.578 


Table 4.4 showing normalized RMS errors of the reconstructed images for elliptical 
Aluminium insert implanted in steel from 4-view projection data with different relaxation 
parameters. 
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Serial Number 

Relaxation 

Number of 

Normalized RMS 


Parameter (A.) 

iterations 

error 

1 

Oo 

439 

18.178 

2 


427 

18.341 

3 


413 

18.564 

4 


395 

18.265 

5 

0.250 

378 

18.127 

6 

0.280 

356 


7 

0.300 

334 


8 

0.345 

313 

18.129 

9 


292 

17.621 

10 

0.350 

261 


11 

0.371 

223 


12 

0.380 

212 

18.674 

13 

0.400 

194 

18.546 

14 


245 

18.451 

15 


349 

18.821 

16 

1.500 

312 

19.869 

17 

1.800 

181 

18.960 


Table 4.5 showing normalized RMS errors of the reconstructed images for square Nickel 
insert implanted in steel from 4-view projection data with different relaxation parameters. 


Serial Number 

Relaxation 
Parameter (A.) 

Number of 
iterations 

Normalized RMS 

error 

1 

0.100 


17.673 

2 


400 

17.341 

3 


385 

17.762 

4 

hhi^sbhhi 

331 

16.897 

5 



17.756 

6 


267 

16.367 

7 

0.400 

224 

17.031 

8 

0.450 

178 

16.987 

9 

0.500 

161 

17.413 

10 

0.550 

149 ^ 

18.098 

11 

0.600 

124 

17.786 

12 

0.650 

112 

18.005 

13 

0.700 

107 

18.737 


Table 4.6 showing normalized RMS errors of the reconstructed images for circular 
Nickel insert implanted in steel from 4-view projection data with different relaxation 
parameters. 
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Serial Number 

Relaxation 
Parameter (X) 

Number of 
iterations 

Normalized RMS 
error ’ 

1 

0.100 

421 

17.112 

2 

0.150 

411 

17.542 

3 


389 

17.983 

4 


356 


5 

0.300 

311 


6 

0.350 

278 


7 

0.400 

233 

16.945 

8 


209 

17.114 

9 

0.500 

187 

17.453 


0.550 

168 

18.213 

11 

r 0.600 

145 

17.981 

12 

0.650 

123 

17.969 

13 

0.700 

119 

18.236 


Table 4.7 showing normalized RMS errors of the reconstructed images for elliptical 
Nickel insert implanted in steel from 4-view projection data with different relaxation 
parameters. 


Serial Number 

Relaxation 
Parameter (X) 

Number of 
iterations 

Normalized RMS 

error 

1 

0.100 

476 

18.213 

2 

0.150 

453 

18.986 

3 

0.190 

411 

18.453 

4 

0.200 

397 

18.342 

5 

0.250 

376 

19.112 

6 

0.280 

361 

18.897 

7 

0.300 

316 

17.985 

8 

0.345 

297 

18.156 

9 

0.348 

286 

18.211 

10 

0.350 

265 

17.967 

11 

0.371 

243 

18.673 

12 

0.380 

221 

19.096 

13 

0.400 

197 

18.543 

14 

0.500 

231 

19.003 

15 

1.000 

344 

19.311 

16 

1.500 

326 

19.278 

17 

1.800 

201 I 

18.891 


Table 4.8 showing normalized RMS errors of the recoristructed images for square Brass 
insert implanted in steel from 4-view projection data with different relaxation parameters. 
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Serial Number 

Relaxation 
Parameter (A,) 

Number of 
iterations 

Normalized RMS 
error i 

1 

0.100 

432 

17.245 

2 

0.150 

419 

17.678 

3 


396 


4 

0.250 

371 

16.983 

5 

0.300 

341 


6 




7 


276 

17.112 

8 


245 

17.456 

9 

0.500 

207 

17.765 

10 

0.550 

189 

18.115 

11 

0.600 

156 

18.167 

12 

0.650 

131 


13 

0.700 

120 

18.184 


Table 4.9 showing normalized RMS errors of the reconstructed images for circular Brass 
insert implanted in steel from 4-view projection data with different relaxation parameters. 


Serial Number 

Relaxation 
Parameter (A) 

Number of 
iterations 

Normalized RMS 
error 

1 

0.100 

443 

17.111 

2 

0.150 

431 

17.786 

3 

0.200 

404 

17.984 

4 

0.250 

398 

17.009 

5 

0.300 

367 

17.650 

6 

0.350 

321 

17.116 

7 

0.400 

287 

16.981 

8 

0.450 

273 

17.237 

9 

0.500 n 

231 

17.874 

10 

0.550 

212 

18.278 

11 

0.600 

189 

18.347 

12 

0.650 

157 

17.783 

13 

0.700 

134 

18.229 


Table 4.10 showing normalized RMS errors of the reconstructed images for elliptical 
Brass insert implanted in steel from 4-view projection data with different relaxation 
parameters. 
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Serial Number 

Relaxation 
Parameter (X) 

Number of 
iterations 

Normalized RMS 
error ' 

1 

0.100 

443 

21.211 

2 

0.150 

431 

20.734 

3 

. 0.200 

404 

19.453 

4 

0.250 

398 

19.006 

5 

0.300 

367 

20.451 

6 

0.350 

321 

21.115 

7 

0.400 

287 

21.145 

8 

0.450 

273 

20.931 

9 

0.500 

231 

19.567 

10 

0.550 

212 

19.341 

11 

0.600 

189 

19.113 

12 

0.650 

157 

21.340 

13 

0.700 

134 

21.098 


Table 4.11 showing normalized RMS errors of the reconstructed images for circular 
Brass insert implanted in steel from 4-view projection data with different relaxation 
parameters using straight ray assumption. The normalized errors are consistently higher 
in this case suggesting poorer reconstruction 



Serial Number 

Relaxation 
Parameter (1) 

Number of 
iterations 

Normalized RMS 

error 

1 

0.100 

411 

20.156 

2 

0.150 

402 

19.113 

3 

0.200 

389 

21.453 

4 

0.250 

378 

21.005 

5 

0.300 

361 

20.568 

6 

0.350 

343 

19.456 

7 

0.400 

292 

21.317 

8 

0.450 

265 

20.561 

9 

0.500 

223 

19.908 

10 

0.550 

207 

20.231 

11 

0.600 

191 

21.341 

12 

0.650 

151 

20.567 

13 

0.700 

137 

20.213 


Table 4.12 showing normalized RMS errors of the reconstructed images for elliptical 
Brass insert implanted in steel from 4-view projection data vwth different relaxation 
parameters using straight ray assumption. The normalized errors are consistently higher 
in this case suggesting poorer reconstruction 


Number of iterations 

Running time for the code 
(straight ray approach) 

Running time for the code 
(curved ray approach) 

50 

20 minutes 

- 2 hours 

100 

45 minutes 

4 hours 

150 

60 minutes 

5.5 hours 


90 minutes 

8 hours 

250 

120 minutes 

10 hours 


Table 4.13 showing a comparison of the times taken to execute the code with and without 
ray bending for various iterations 









Figure 4.3 Ray-paths for three 90° fan-beams, one passing through the insert, the other 
two just lying on either side of the insert. The plot is for circular Aluminium insert 



Figure 4.4 Ray-paths for three 0° fan-beams, one passing through the insert, the other two 
just lying on either side of the insert. The plot is for circular Nickel insert 
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Figure 4.5 Ray-paths for three 45° fan-beams, one passing through the insert, the other 
two just lying on either side of the insert. The plot is taken for circular Brass insert after 5 
iterations. 



Figure 4.6 Typical ray-paths for circular Aluminium insert implanted in Steel for 
relaxation parameter = 0.1. The plot obtained is after 10 iterations. 
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Figure 4.7 Typical ray-paths for circular Nickel insert implanted in Steel for relaxation 
parameter = 0.1. The plot is after 10 iterations. 
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Figure 4.8 Plot of normalised RMS errors vs number of iterations with and \vithout ray 

bending for circular Brass insert 



Figure 4.9 Plot of normalised RMS errors vs number of iterations with and without ray 

bending for elliptical Brass insert 
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Some images are shown below for the case of Brass insert in Steel where straight ray 
assumption has been taken. The images are shown for the purpose of comparison with 
images obtained using curved ray approach. 



Figure 4.1 0 Tomographic reconstruction for square Brass 
insert impianted in Steel with straight ray assumption for 
reiaxation parameter = 0.1 



Figure 4.11 Tomographic reconstruction for circular Brass 
insert implanted in Steel with straight ray assumption for 
relaxation parameter = 0.1 



Figure 4.12 Tomographic reconstruction for elliptical Brass 
insert with straight ray assumption for relaxation parameter = 0.1 



The reconstructed images for the square shaped Aluminium insert implanted in Steel for 
four different relaxation parameters are given in figure 4. 1 3 *-4. 1 6. ' 



Figure 4.1 3 Tomographic reconstruction for square Aluminium Figure 4.1 4 Tomographic reconstruction for square Aluminium 
insert implanted in Steel for relaxation parameter = 0.1 insert implanted in Steel for relaxation parameter = 0.2 



Figure 4.15 Tomographic reconstruction for square Aluminium 
insert implanted in Steel for relaxation parameter = 0.3 


Figure 4.16 Tomographic reconstruction for square Aluminium 
insert implanted in Steel for relaxation parameter = 0.4 



Figure 4.1 9 Tomographic reconstruction for circular Aluminium Figure 4,20 Tomographic reconstri^ton for circular Aluminium 
insert implanted In Steel for relaxation parameter = 0.3 insert implanted in Steel for relaxation parameter = 0.4 
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Figure 4.23 Tomographic reconstruction for elliptical Aluminium Figure 4.24 Tomographic reconstructton for eiptical Aluminium 

insert implanted in Steel for relaxation parameter = 0.3 insert implanted in Steel for relaxation parameter = 0.4 

Wife So A 
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Figure 4.25 Contour plot for square Aluminium insert in Steel 
The values in the l^end are the normalized slowness values 
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Figure 4.26 Contour plot for circular Aluminium insert in Steel 
The values in the legend are the normalized slowness values 






0,84 

0.87 

0^ 

0,90 


Figure 4,27 Contour plot for elliptical Alumniurn insert in Steel 
The values in the legend are the normalized slowness values 




The reconstructed images for the square shaped Nickel insert implanted in Steel for four 
different relaxation parameters are given in figure 4.28 -4.31. 



Figure 4.28 Tomographic reconstruction for square Nickel 
insert implanted in Steel for relaxation parameter = 0.1 



Figure 4.30 Tomographic reconstruction for square Nickel 
insert implanted in Steel for relaxation parameter = 0.3 



Figure 429 Tomographic reconstruction for square Nickel 
insert implanted in Steel for relaxation parameter = 02 



Figure 4.31 Tomographic reconstruction for square Nickel 
insert implanted in Steel for relaxation parameter = 0.4 
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The reconstructed images for the circular shaped Nickel insert implanted in Steel for four 
different relaxation parameters are given in figure 4.32 - 4.35. 



Figure 4.32 Tomographic reconstruction for circular Nickel Figure 4.33 Tomographic reconstruction for circular Nickel 

insert implanted in Steel for relaxation parameter = 0.1 insert implanted in Steel for relaxation parameter = 0.2 




Figure 4.34 Tomographic reconstruction for circular Nickel 
insert implanted in Steel for relaxation parameter = 0.3 



Figure 4.35 Tomographic reconstruction for circular Nickel 
insert implanted in Steel for relaxation parameter = 0.4 
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The reconstructed images for the elliptical shaped Nickel insert implanted in Steel for 
four different relaxation parameters are given in figure 4.36 - 4.39. 



Figure 4.38 Tomographic reconstruction for elliptical Nickel 
insert implanted in Steel for relaxation parameter = 0.3 


Figure 4.39 Tomographic reconstruction for elliptical Nickel 
insert implanted in Steel for relaxation parameter = 0.4 












Figure 4.42 Tomographic reconstruction for square Brass 
insert implanted in Steel for relaxation parameter = 0.3 


Figure 4.43 tomographic reconstruction for square Brass 
insert implanted in Steel for relaxation parameter = 0.4 






The reconstructed images for the elliptical shaped Brass insert implanted in Steel for four 
different relaxation parameters are given in figure 4.48 -4.51. 



Figure 4.48 Tomographic reconstruction for elliptical Brass 
insert implanted in Steel for relaxation parameter = 0.1 



Figure 4.50 Tomographic reconstruction for elliptical Brass 
insert implanted In Steel for relaxation parameter = 0.3 



Figure 4.49 Tomographic reconstruction for elliptical Brass 
insert implanted in Steel for relaxation parameter = Q2 



Figure 4.51 Tomographic reconstruction for elliptical Brass 
insert implanted in Steel for relaxation parameter = 0.4 
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4.6 Closure 


Reconstructed images have been shown for all the three combinations of parent material 
and insert material. The quality of reconstruction is quite satisfactory for Steel- 
Aluminium pair and also for Steel-Nickel pair. However, it is not as good for Steel-Brass 
pair because of the larger velocity difference between Steel and Brass. The quality of 
reconstruction is also foxmd to be best for relaxation parameters varying between 0.1 and 
0.4. Hence all the images are shown for relaxation parameters lying within this range. 
Typical ray paths have been shown after the completion of a few iterations. Finally, a 
number of tables have been provided which furnish the variation in the RMS error as the 
relaxation parameter is varied. A table showing the comparison of the times taken for the 
algorithm to converge with and without ray bending has also been provided. 
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Chapter 5 


Conclusion and scope of future work 

5.1 Conclusion 

In the present work, standard techniques have been applied for the tomographic 

reconstruction of samples with artificially implanted insertions. Based on the results 

obtained, the following inferences can be drawn: 

• For isotropic materials the acoustic analogue of Snell’s law is used to predict the path 
of the ultrasonic ray with reasonable accuracy in the presence of a refraction artifact. 

• For tomographic reconstruction of circular and elliptical defects, recourse has been 
taken of standard edge detection techniques to find out correctly the interface for 
refraction. Without edge detection, the errors sometimes become unacceptably high. 

• The normalized RMS error, which is a measure of the accuracy of reconstruction, is 
found to be consistently higher for images which have been reconstructed without 
using curved ray approach. The normalized RMS errors lie in the range of 21% for 
images obtained using straight ray approach and 18% for images obtained using 
curved ray approach. 
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